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Abstract 

We consider the problem of optical tomographic imaging in the mesoscopic regime where the 
photon mean free path is of order of the system size. Within the accuracy of the single-scattering 
approximation to the radiative transport equation, we show that it is possible to recover the extinc- 
tion coefficient of an inhomogeneous medium from angularly-resolved measurements. Applications 
to biomedical imaging are described and illustrated with numerical simulations. 



1 



I. INTRODUCTION 



There has been considerable recent interest in the development of experimental meth- 
ods for three-dimensional optical imaging of biological systems. Applications range from 
imaging of optically-thin (by which we mean nearly transparent) cellular and sub-cellular 
structures to optically-thick systems at the whole organ level in which multiple scattering 
of light occurs. In optically-thin systems, confocal microscopy [l| and optical coherence 
otoscopy 3 can be used to generate three-dimensional image, by optical sectioning Al- 
ternatively, computed imaging methods such as opticalprojection tomography (OPT) |3[ or 
interferometric synthetic aperture microscopy (ISAM) [4], Q may be employed to reconstruct 
three-dimensional images by inversion of a suitable integral equation. In the case of OPT, 
the effects of scattering are ignored and geometrical optics is used to describe the propa- 
gation of light. The sample structure is then recovered by inversion of a Radon transform 
which relates the extinction coefficient of the sample to the measured intensity of the optical 
field. In the case of ISAM, the effects of scattering are accounted for within the accuracy of 
the first Born approximation to the wave equation. An inverse scattering problem is then 
solved to recover the susceptibility of the sample from interferometric measurements of the 
cross-correlation function of the optical field. 

In optically-thick systems, multiple scattering of the illuminating field creates a funda- 
mental obstruction to image formation. If the medium is macroscopically large and weakly 
absorbing, only diffuse light is transmitted. By making use of the diffusion approximation 
(DA) to the radiative transport equation (RTE) and solving an appropriate inverse problem, 
the aforementioned difficulty may, to some extent, be overcome. This approach forms the 
basis of diffusion tomography which can be used to reconstruct images with sub- centimeter 
resolution of highly-scattering media such as the human breast 6|. The relatively low quality 
of reconstructed images is due to severe ill-posedness of the inverse problem. 

Despite significant recent progress in optical imaging of both optically thin and thick 
media, little has been done for imaging of systems of intermediate optical thickness. This 
represents the subject of this study. In radiative transport theory, such systems are referred 
to as mesoscopic, meaning that the photon mean free path (also known as the scattering 
length) is of order of the system size |7| ■ In the mesoscopic scattering regime, applications to 
biological systems include engineered tissues and semitransparent organisms such as zebra 
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fish, animal embryos, or small-animal extremities. In this case, light exhibits sufficiently 
strong scattering so that the image reconstruction methods of computed tomography are 
not applicable, yet the detected light is not diffuse and the diffusion tomography can not be 
employed. 

In mesoscopic systems, the DA does not hold and the RTE must be used to describe the 
propagation of light In this study, light transport in the mesoscopic regime is described 
by the first-order scattering approximation to the radiative transport equation. This en- 
ables the derivation of a relationship between the extinction coefficient of the medium and 
the single-scattered light intensity, which represents the basis for a novel three-dimensional 
optical imagining technique that we propose and refer to as single scattering optical tomog- 
raphy (SSOT). SSOT uses angularly-selective measurements of scattered light intensity to 
reconstruct the optical properties of macroscopically inhomogeneous media, assuming that 
the measured light is predominantly single-scattered. The image reconstruction problem 
of SSOT consists of inverting a generalization of the Radon transform in which the inte- 
gral of the extinction coefficient along a broken ray (which corresponds to the path of a 
single-scattered photon) is related to the measured intensity. 

Our results are remarkable is several regards. First, similar to the case of computed 
tomography, inversion of the broken-ray Radon transform is only mildly ill-posed. Second, 
the inverse problem of the SSOT is two-dimensional and three-dimensional image recon- 
struction can be performed slice- by-slice. Third, in contrast to computed tomography, the 
experimental implementation of SSOT does not require rotating the imaging device around 
the sample to acquire data from multiple projections. Therefore, SSOT can be used in the 
backscattering geometry. Finally, SSOT makes use of intensity measurements, as distinct 
from the more technically challenging experiments of optical coherence microscopy or ISAM, 
which require information about the optical phase. 

This paper is organized as follows. In Sec. HIl we introduce the single-scattering ap- 
proximation appropriate for the mesoscopic regime of radiative transport. We then derive 
a relationship between the scattering and absorption coefficients and the single-scattered 
intensity. This relationship is then exploited in Sec. II III to discuss the physical principles 
of SSOT. In Sec. IIV} numerical algorithms for both the forward and inverse problems are 
presented and illustrated in computer simulations. 
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II. MESOSCOPIC RADIATIVE TRANSPORT 



We begin by considering the propagation of light in a random medium of volume V. The 
specific intensity I(r, s) is the intensity measured at the point r and in the direction s, and 
is assumed to obey the time-independent RTE 

[s • V + Mr) + /u s (r)] J(r, s) = /i s (r) J A(a, s')/(r, s')d 2 s' , r e V , (1) 

where /x a ( r ) an d /-* s ( r ) are the absorption and scattering coefficients. The phase function 
A(s, s') describes the conditional probability that a photon traveling in the direction s is 
scattered into the direction s' and is normalized so that J A(s, s')d 2 s' = 1 for all s. Eq. ([I]) 
is supplemented by a boundary condition of the form 

J(r,s) =/inc(r,s) , s-n(r)<0, r G dV , (2) 

where n is the outward unit normal to dV and / inc is the incident specific intensity at the 
boundary. 

The RTE (CQ) together with the boundary condition (T5]) can be equivalently formulated 
as the integral equation 

/(r, s) = / 6 (r, s) + J G(r, s; r', s')^(r')A(s\ s")/(r, s")d 3 r'd 2 s'd 2 s" . (3) 

Here If, is the ballistic (unscattered) contribution to the specific intensity and G is Green's 
function for the ballistic RTE, which satisfies the equation 

[s- V + // a (r) + /x fl (r)]/ 6 (r,s) = , (4) 

and obeys the boundary condition ([2]). If a narrow collimated beam of intensity Iq is incident 
on the medium at the point r x in the direction s 1; then /^(r, s) is given by 

J 6 (r,s) = I G(r 1 s;r ll s 1 ) , (5) 
where ballistic Green's function G(r, s; r', s') is expressed as 



r — r' 



G(r, s; r', s') = ^(r, r')6 ( s' - -— - { ) 6{s - s'). (6) 
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Here 

1 



3( r > r ') = | r _ r ,|2 eX P 







» t [r' + ly-^)d£ , (7) 

and the Dirac delta function 5(s — s ) is defined by 

5(s — s ) = — </?§' )^ ( cos ^§ — cos 0§' ) , (8) 

where we have introduced the extinction (attenuation) coefficient fit — A*o + fts, and 6* and 
<£> are the polar angles of the respective unit vectors. Note that g is the angularly-averaged 
ballistic Green's function, 



g(r,r') = j d 2 sd 2 s'G(r, s; r', §') . (9) 

By iterating Eq. ([3]) starting from = 1^, corresponding to ballistic light propagation, 
we obtain 

J(r,s)=/W(r,s)+/( 1 )(r,s) + /( 2 )(r, §) + ■•■ , (10) 
where each term of the series is given by 

/W(r, g) = y d 3 r'd 2 s'd 2 s"G(r, s; r', s> s (r')A(s', s")/^ -1 ^ s"). (11) 

The Born series ( jTOl can be regarded as an expansion in the number of scattering events, 
each term corresponding to a successively higher order of scattering [q| . The convergence of 
this series requires that 1 1 J d 3 r'd 2 s'G(r, s; r', s')/j, s (r')A(s', s") | |oo < 1, where 1 1 * j |oo is the L°° 
norm. For an isotropically scattering random medium, this norm can be calculated to be of 

n 

the order [9J (/x s //^)(l — exp(— /x t i2)), where R is the characteristic size of the system. Thus, 
the convergence requirement for the Born series associated with the RTE is always satisfied. 
For the system under investigation, this norm varies from 0.16 to 0.5, as the amount of 
scattering is increased such that n s R varies from 1.6 to 6.4, respectively. Therefore, very 
rapid convergence is expected. 

In general, the specific intensity can be decomposed as / = Ib+/ S , where I s is the scattered 
part of the specific intensity. Within the accuracy of the single-scattering approximation, I s 
is given by the expression 



^(r,s) 



J d 3 r'd 2 s'd 2 s"G(r 1 s; r, s')fi s (r')A(s' , s")I b (r', s") . (12) 
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FIG. 1: (Color online) Illustrating the geometrical quantities used in Eq. (I13|) . BR denotes broken 
ray. 

Note that the single scattering approximation is expected to hold in the mesoscopic regime 
of radiative transport, where the system size is of order the scattering length (defined to be 

We now assume that the sample is a slab of width L and that the beam is incident on 
one face of the slab at the point ri in the direction §1 and that the transmitted intensity 
is measured on the opposite face of the slab at the point r 2 in the direction s 2 (see Fig. [1]). 
We denote the scattered intensity measured in an such experiment by I s (ri, §1; r 2 , s 2 ). Per- 
forming the integral in expression (112j) with lb given by (jSJ) and using ((6]) yields 



J s (ri, §i; r 2 , s 2 ) = IoQ(n -Q\- 9 2 )S (|^ §1 - (p 



S2 



-7T 



r 21 sin 6*i sin 9 2 



x exp 



-L\ />L2 

(13) 



/ /i t (ri + tsi)M - / /ii(R 2 i + £s 2 
Jo Jo 



where Q(x) is the step function, R 2 i is the position of the turning point of the ray, r 2 i = 
r 2 — ri , r 21 = |r 21 |, Li = |R 21 — ri|, L 2 = |r 2 — R 2 i|, and the angles Q\ and 6 2 are defined by 
cos#i i2 = r 2 i ■ Si 2 . The details of the derivation of Eq. ffTB"]) are presented in the Appendix. 
We note the following important relations: 

R21 = ri + = r 2 - L 2 s 2 (14) 

sin(6>i + 6*2) sin(6'i + fc> 2 ) 

The physical meaning of the various terms in (|T3|) is as follows. First, the angles </9g 12 in 
the Dirac delta function 8 (|<£>§i — ^§2! — t 1 ") are the azimuthal angles of the unit vectors Si j2 
in a reference frame whose z-axis intersects both the position of the source and the detector 
(this axis is shown by a dashed line in Fig. [1] and should not be confused with the z-axis 
of the laboratory frame shown by a solid line). The presence of this one-dimensional delta 
function is the manifestation of the fact that two straight rays exiting from the points r x 
and r 2 in the directions §1 and — s 2 , respectively, can intersect only if §1 and s 2 and r 21 are 
in the same plane (equivalently, if tp^ — ifs 2 = 0, ±7r) and point into different half-planes 
(this requires that ip^ — V 9 ^ = ±tt)- Second, the point of intersection exists within the plane 
only if 61+62 < 7r, which is expressed by the step function 0(ir — 6\ — 6* 2 ). We note that if, 
additionally, §1 and s 2 are restricted so that z-§i > and z-s 2 < (si points into the slab and 
s 2 points out of the slab), then R 21 lies within the slab. Third, the factor /i s (R 21 )y4(s 2 , §1) 
is the "probability" that the ray is scattered at the point r = R 2 i and changes direction 
from §1 to s 2 . This factor is, in general, position-dependent. Fourth, l/r 2 i sin#i sin# 2 is 
a geometrical factor. We note that it can be equivalently rewritten as r-ixj H-zxHvii where 
H 2 i and H12 are the two heights of the triangle (r 1; R 2 i, r 2 ) drawn from the vertices ri and 
r 2 , respectively, as shown in Fig. [TJ Finally, it can be seen that the integral of fi t i n the 
argument of the exponential is evaluated along a broken ray which begins at ri, travels in 
the direction §1, turns at the point R 2i , travels in the direction s 2 , and exits the slab at r 2 . 
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III. PHYSICAL PRINCIPLES OF SSOT 

The physical principle of SSOT is illustrated in Fig. [2] where a slab-shaped sample is 
illuminated by a normally incident beam. In the absence of scattering, the beam would 
propagate ballistically, as shown by the green ray. Detection of such unscattered rays is the 
basis of computed tomography. In the presence of scattering, the ray can change direction 
as shown in Fig. |2(a). Of course, scattering does not result in elimination of the ballistic 
ray. However, it is possible to employ angularly- selective detectors that do not register 
the ballistic component of the transmitted light. In the example shown in Fig. [2](a), it is 
assumed that only the intensity of the broken ray shown by the red line is detected. To 
avoid detection of the ballistic component of the transmitted light, the angularly-selective 
source and detector are not aligned with each other. Moreover, the data can be collected 
either on opposite sides of the slab (transmission measurements), or in the backscattering 
geometry. In both cases, rotation of the instrument around the sample is not required. 

By utilizing multiple incident beams and detecting the light exiting the sample at different 
points, as shown in Fig. [2]^b), we will see that it is possible to collect sufficient data to 
reconstruct the spatial distribution of the attenuation coefficient in a fixed transverse slice 
of the slab. In addition to varying the source and detector positions, one can also vary the 
angles of incidence and detection. In principle, this can provide additional information for 
simultaneous reconstruction of absorption and scattering coefficients of the sample, a topic 
we will consider elsewhere. 

The image reconstruction problem of SSOT is to reconstruct fi a , fi s and A from measure- 
ments of I s . For simplicity, we assume that /i s and A are known, in which case we wish to 
determine /i a . To proceed, we use Eq. (1131) to separate the known or measured quantities 
from those that need to be reconstructed, and introduce the data function 



where we have assumed that A is position-independent. Note that if I s is experimentally 
measured, the angular integration on the right-hand side of ( fl6l) does not need to be per- 
formed numerically. The measured data is necessarily integrated in a narrow interval of 
(fis 2 due to the finite aperture and acceptance angle of the detector. Note also that the 
above definition is only applicable for such configurations of sources and detectors such that 




(16) 
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FIG. 2: (Color online) (a) Schematic illustration of the broken-ray (or single-scattered ray, denoted 
here by SSR) trajectory, (b) Source-detector arrangement for SSOT. Reconstruction is performed 
independently slice-by slice. The blue rectangle represents the area in which a reconstruction can 
be performed. 

61+62 < tv. Otherwise, any measured intensity is due to higher-order terms in the collision 
expansion which are not accounted for ( fl~3l) . 

Making use of the definition ( |T6l) and Eq. ( jl~3l) . we find that \i t obeys the integral equation 

0(r 1 ,s 1 ;r 2 ,s 2 )= / ■ (17) 

J BR(ri,si;r2,S2) 

Here the integral is evaluated along the broken ray (such as the one shown in Fig. []]), which 
is uniquely denned by the source and detector positions and orientations, and I is the linear 
coordinate along the ray. 

According to (|T7|) . the attenuation function is linearly related to the data function. In 
this respect, the mathematical structure of SSOT is similar to the problem of inverting the 
Radon transform in computed tomography except that the integrals are evaluated along 
broken rays. Since \xt may be regarded as a function of two variables, it is sufficient to 
consider only two-dimensional measurements. One possible choice is to vary the source and 
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detector coordinates, y\ and y 2 , while keeping the angles of incidence and detection fa and 
fa fixed. By fa and fa we mean here the angles between the z-axis of the laboratory frame 
and the unit vectors §1 and s 2 , respectively. Note that these angles are not equal to the 
angles 9i and 82 shown in Fig. [TJ The latter can vary in the measurement scheme described 
in this section, while fa and fa are fixed. Below, we omit fa and fa from the list of formal 
arguments of the data function and consider the equation 

0(2/1,2/2)= / »t(y{£),z{e))de, (18) 

JBR{y im ) 

where 0(2/1,2/2) is the two-dimensional data function. 

As explained above, the selection of the points and directions of incidence and detection 
define a slice in which fi t is to be reconstructed. In Fig. [TJ this slice coincides with the yz- 
plane of the laboratory frame. Assuming that the ^-coordinate is fixed, we then regard fit 
as a function of y and z. Three-dimensional reconstruction is then performed slice-by-slice. 



IV. IMAGE RECONSTRUCTION 



In this section we illustrate image reconstruction in SSOT using a numerical technique 
based on discretization and algebraic inversion of the two-dimensional integral equation 
( TIB]) . We note that more sophisticated image reconstruction procedures which utilize the 
translational invariance of rays may also be derived. These methods are conceptually similar 
to those previously developed for optical diffusion tomography 
elsewhere. 
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ll| and will be described 



A. Forward Problem 

We begin by describing a method to generate simulated forward data to test the SSOT 
image reconstruction. Assuming an isotropically scattering sample with A(s, s') = 1/47T, it 



can be shown from ([3D 
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13] that the specific intensity everywhere inside the sample is 



related to the density of electromagnetic energy u(r) = J I(r, s)d 2 s by the formula 



J(r, s) = I b (r, s) + 



— I G(r,s;r / ,s / )w a (r')«(rOd 3 r / dV, 
4tt J 



(19) 



where u(r) satisfies the integral equation 
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u(r) = u b (r) + 



1 



g(r,r')fi s (r')u(r')d 3 r' , 



(20) 



Here g(r, r') is given by ([7]) and u b (r) = f I b (r, s)d 2 s is the ballistic energy density We note 
that the assumption of isotropic scattering is the most stringent test for SSOT. 

The specific intensity is computed by first solving the integral equation (1201) and then 
substituting the obtained solution u(r) into ( |T9l) (where the ballistic part I b may be ignored). 
Note that I(r, s) calculated from (Tl9|) satisfies the boundary conditions at all surfaces. We 
also stress that this numerical approach is non-perturbative and includes all orders of scat- 
tering. 

In the simulations shown below, ( !20l is discretized on a rectangular grid and solved by 
standard methods of linear algebra. The energy density u(r) is assumed constant within 
each cubic cell, and the corresponding values u n = w(r n ), where r n is the center of the nth 
cubic cell, obey the algebraic system of equations obtained by discretizing Eq. (1201 . The 
off-diagonal elements of the matrix of this system corresponding to the integral on the right- 
hand side of (1201 are given by (/i s /i 3 /47r)g(r m , r n ), where h is the discretization step. Here we 
can take advantage of the fact that fi s was set to be constant throughout the sample, while 
the inhomogeneities were assumed to be purely absorbing. Computation of the diagonal 
elements is slightly more involved because g(r, r') diverges when r — > r'. In this case, we 
need to find an approximation for the integral 



where the integration is carried out over the nth cell. While integration over a cubic volume 
is difficult, the important fact is that the singularity in g(r, r') is integrable. We then write, 
approximately, 



where R eq = (3/47r) 1 / 3 /i is the radius of a sphere of equivalent volume. For a sufficiently fine 
disctretization, /it-R C q *C 1, which allows us to write g(0, r) w 1/r 2 . This leads to S = fi s R e q, 
and the discretized version of Eq. (1201 becomes 




(21) 




(22) 



(1 - n s R eq )u, 



4tt 



^2 #( r r» r m )u m = u b {r n ) . 



(23) 



11 



We note that this set of equations is an accurate approximation to the integral equation 
(120]) if fi s R eq ~ /i s /i <C 1. However, since in practice this inequality may be not very strong, 
the term fi s R eq on the left-hand side of (12UI) is not neglected. 
Eq. (T2"5j) can be written in matrix notation as 



W\u) = \b) (24) 



where 



(n\W\m) = -5 nm 
a 



h 2 



|2 



exp 



(M (r m + u(r 



(1 - 6 nm ) , (25) 



with 5 nm being the Kronecker delta function, a being a dimensionless coupling constant 
defined by 

" = 4,(1 - • (26 » 

(n|w) = w n , and (n\b) = 4TTUh(r n ) / fi s h. Note that the quantity «&(r n ) is defined as an average 
over the nth cell, namely, Ub(r n ) = h~ 3 j v Ub(r)d 3 r. 

Computing the elements of the matrix W requires the evaluation of the integrals on the 
right-hand side of (1251) . For a homogeneous sample, this can be performed analytically. For 
an inhomogeneous sample, it is done numerically. Once this is accomplished, (1241) is then 
solved by matrix inversion. We note that the symmetric matrix W is well-conditioned fl^ . 
This is illustrated in Fig. [3] where we plot all the eigenvalues of W for /x s = 0.08/i _1 with 
the set of absorbing inhomogeneities described in section IIVCI We also note that W is 
positive-definite, so that u and u& are positive. Also, it was verified in all simulations that 
the diffuse component of the density, defined as the quantity u — Ub, was positive everywhere 
inside the sample. 

Special attention must be given to the effects of discretization when computing the single- 
scattered intensity, the data function, and the forward data. Since the computations in- 
volve discrete rays, employing the expressions (fT3|) and (fl6l) that contain the delta-function 
S (!<£>§! —V 9 ^ I — tt) an d the geometrical factor I/V21 sin 82 sin 9\ becomes cumbersome. In- 
stead, we have derived the discrete analogues of these expressions starting form Eq. (1121) 
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FIG. 3: (Color online) All 53, 680 eigenvalues w n of the matrix W of Eq. [2H computed for fj, s = 
0.08/1" 1 and the set of absorbing inhomogeneities described in section HV Bl 

from which (|T3|) has been obtained. In particular, the expression for the single-scattered 
intensity (the discrete analog of ffT3"|) ) is obtained from Eq. ffT2"j) as 

a s h 3 

I s (r 2 , s 2 ; n, §i) = -J— g(r 2 , R 2 iK(R 2 i) , (27) 

47T 

where Ub(R. 2 i) is the average of u&(r) over the cell that contains R 2 i. This expression and 
the expression (j5J) for the ballistic intensity suggest the definition of the data function of the 
form 



(r 2 ,s 2 ;ri,si) = -In 



47r/ s (r2,s 2 ;ri,s/ 



(28) 



h 3 I Q fi a 

This is the discrete analogue of (fl~6j) . Finally, the data function is calculated according with 
this expression and using the specific intensity obtained from the discretized version of (|T9l) . 

J(r 2 ,s 2 ;ri,Si) = -|— g(r 2 ,r n )u n , (29) 

r2"r„=s 2 |r 2 -r ra | 

where u n must be computed numerically for the selected source. The condition on the 
sum means that summation is performed only over cells that are intersected by the ray 
exiting from the detection point r 2 in the direction s 2 . The above formula is valid for the 
specific measurement scheme which obtains when the intersection length of all such rays with 
any cubic cell is constant. Otherwise, a more complicated numerical integration must be 
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employed. We note while Eqs. (125]) and ([29]) are employed for simulated data, Eqs. ([TBI and 
([16]) must be used when experimental data is available. Also, the use of experimental data 
avoids many mathematical complications that arise due to discretization of rays, needed to 
solve the forward RTE problem, usually computationally intensive. 

In order to model the noise in measured data, the specific intensity obtained from the 
forward solver was scaled and rounded so that it was represented by 16-bit unsigned integers, 
similar to the measurement by a typical CCD camera. Then a statistically-independent 
positively-defined random variable was added to each measurement. The random variables 
were evenly distributed in the interval [0,n/ av ], where n is the noise level indicated in the 
figure legends below and I m is the average measured intensity (a 16-bit integer). The DC part 
(the positive background) of the intensity was not subtracted (this procedure is commonly 
applied to the digital output of CCD chips). Then the simulated intensity measurements, 
together with the appropriately scaled incident intensity I were substituted into ([16]) to 
obtain the data function <fi. 

B. Inverse Problem 

We now describe the method by which we invert the integral equation (TT8]) . The discrete 
version of ( 118]) has the form 



where the same grid is employed as is used to solve the forward problem, £ vn is the length 
of the intersection of the broken ray indexed by v = (yi, y 2 ) with the nth cubic cell (located 
within the selected x-slice of the sample). The matrix elements £ are determined from 
simple geometric considerations. The matrix form of ( |30]) is 




(30) 



n 



C\lh) = 10) • 



(31) 



The equation (131]) can be solved using a regularized pseudoinverse [15|, namely, 



(32) 



Here (£*£) 1 is understood in the following sense: 



14 



n 



0.01 







200 



400 



600 



800 



1000 



1200 



FIG. 4: (Color online) All 1, 156 non-zero singular values a n of the matrix L denned by Eq. ([30 
(the size of C in this example is 1, 600 x 1, 156). 



|/n>(/n| 



(33) 



where Q(x) is the step function, e is a small regularization parameter, and |/ n ) and a n are 
the singular vectors and singular values [jjj], respectively, of the matrix C. These quantities 
are the solution of the symmetric eigenproblem C*C\f n ) = a^\f n ). A typical spectrum 
of singular values of L for \x s = 0.08/i _1 , 1,600 measurements and 34 2 = 1,156 unknown 
values of [i t (the size of C in this example is 1, 600 x 1, 156, so that the problem is slightly 
overdetermined) is shown in Fig. HI It can be seen that matrix condition number 14j for 
the inverse problem is much larger than for the forward problem. In the example shown in 
Fig. HI the condition number is ~ 10 3 . Thus the inverse problem is very mildly ill-posed. 



C. Numerical Results 

In what follows, we illustrate applications of SSOT to biological imaging. In particular, 
the numerical simulations presented here are relevant for "semi-transparent" systems, such 
as zebra fish or engineered tissues. Also, the physical situations analyzed here are exper- 
imentally encountered for organ tissues at certain wavelengths of the illuminating beam 
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Reconstructions were carried out for a rectangular isotropically scattering sample of di- 
mensions L x = 25h, L y = 122/i and L z = 40h. The background absorption coefficient of 
the sample was equal to 0.01/i _1 and was spatially modulated by absorbing inhomogeneities 
(the target). The target was a set of absorbing inclusions formed in the shape of letters, 
with absorption varying from 0.06/i _1 to 0.2/i" 1 . The inclusions were concentrated only in 
three layers: x = 6h, x = 13h and x = 20h, as shown in the columns marked model of 
Figs. [5J13 The scattering coefficient was constant throughout the sample, with three differ- 
ent values used throughout the simulations corresponding to /i s = 0.04/i _1 , fi s = 0.08/i _1 
and fi s = 0.16/i -1 . Thus, for example, in the case fi s = 0.04/i _1 the contrast of fit (the ratio 
of Ht in the target to the background value) varied from 2.0 in the letters RADIOL to 4.8 
in the letters DEPT. In the case fi s = 0.16/i _1 , the contrast was smaller and varied from 
1.18 to 2.12. We note that the contrast in the total attenuation coefficient depends weakly 
on the background absorption coefficient, compared to its dependence on the background 
scattering coefficient. 

The sources were normally incident on the surface z = 0. The detectors were placed 
on the opposite side of the sample and the specific intensity exiting the surface z = L z at 
the angle of ir/4 with respect to the z-axis was measured. In this situation there are two 
possibilities — the exiting rays either make an angle of 7r/4 or 37r/4 with the y-axis. In some 
cases, data from both directions were used. Note that the distance L z corresponds to the 
slab thickness L used in Sec. [TTJ The optical depth of the sample, fi s L z , varied from 1.6, 
for n s = 0.04/i _1 , to 6.4, for fi s = 0.16/i" 1 . This corresponds to the mesoscopic scattering 
regime in which the image reconstruction method of SSOT is applicable. 

Reconstruction of the total attenuation coefficient fi t was performed in slices x = x s iicc 
separated by a distance Ax = h. For each slice, the source positions were x = x s u cc , y = nh, 
z = 0, with n being integers. The reconstruction area inside each slice was AAh < y < 77h, 
Ah < z < 37h, with the field of view 34/i x 34h. At the noise levels n = and n = 1%, only 
the rays making an angle of 7r/4 with the y-axis were used; for the noise level n = 3%, the 
exiting rays which make an angle of 3-7r/4 with the ?/-axis were also used in order to improve 
image quality of the reconstructions. Also, the regularization parameter e in the regularized 
pseudoinverse (|33|) was varied to obtain the best visual appearance of images. Note that the 
absolute values of the reconstructed fit are not sensitive to the choice of e. Qualitatively, the 
same results are obtained by setting e = 0, although we have found that selecting a small 
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but nonzero value of e tends to slightly improve image quality. 

The results of image reconstructions for various noise levels are presented in Figs. [5]- 
where we show both slices containing inhomogeneities and neighboring slices in which 
inhomogeneities are not present. It can be seen that the spatial resolution depends on the 
noise level and contrast, and can be as good as one discretization step, h. Note that the 
reconstructed images are in very good quantitative agreement with the model (all panels in 
each figure are plotted using the same color scale) and stable in the presence of noise. When 
H s = 0.16/i _1 (Fig. [JJ the optical depth of the sample is fi s L z = 6.4. This is a borderline case 
when scattering is sufficiently strong so that the single-scattering approximation of SSOT 
may be expected to be inaccurate. Indeed, the image quality in Fig. [JJ is markedly worse 
than in Figs. and El yet the letters in the image remain legible. 

We emphasize that the reconstructed images presented here are based on simulated data 
obtained by solving the RTE exactly, thereby accounting for all orders of scattering. For 
samples which are optically thick, the resulting reconstructions evidently exhibit artifacts 
due to the breakdown of the single scattering approximation (which is not possible physically 
but achievable in simulations). Such would be absent if only single-scattered light were 
detected. To illustrate this idea, we present in Fig. [8] reconstructed images for the case 
fi s = 0.16/i^ 1 using forward data in which only single-scattered light is retained. Here, 
instead of solving the full RTE, the data function was directly calculated from ffTTl) . The 
aforementioned procedure of generating data overestimates the performances of the imaging 
method. It can be seen from Fig. 8 that, as expected, significantly better reconstructions 
result. 

V. DISCUSSION 

We have investigated the problem of optical tomography in the mesoscopic regime. 
Within the accuracy of the single scattering approximation to the RTE, we have derived a 
relation between the absorption and scattering coefficients and the specific intensity. In par- 
ticular, for a homogeneously scattering medium, we have shown that the intensity measured 
by an angularly-selective detector is related to the integral of the attenuation coefficient along 
a broken ray. By inverting this relation, we are able to recover the attenuation coefficient of 
the medium. 
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FIG. 5: (Color online) Image reconstruction for a slab with fi s = 0.04/i _1 and various noise levels 
n. The absorbing inhomogeneities are placed in the slices x = 6h, x = 13h and x = 20h, and the 
rows show the slices x = 6h, 7h, 12h, 13/i, 14/i, 19h, and 20h. The same color scale is used for 
all slices, with the maximum (white) corresponding to [it = 0.24/1" 1 and the minimum (black) to 
fit = 0. Red regions correspond to negative values of the reconstructed extinction coefficient. 
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FIG. 6: (Color online) Image reconstruction for a slab with fi s = 0.08/1" 1 . The same color scale 
is used for all slices, with the maximum (white) corresponding to fit = 0.28/i _1 and the minimum 
(black) to fit = 0. All the other details are as for Fig. [5j 
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FIG. 7: (Color online) Image reconstruction for a slab with fi s = 0.16/1" 1 . The same color scale 
is used for all slices, with the maximum (white) corresponding to fit = 0.36/i _1 and the minimum 
(black) to fit = 0. All the other details are as for Fig. [5j 
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FIG. 8: (Color online) Image reconstruction for fi s = O.I6/1 , for a data function corresponding 
only to single-scattered light and calculated according to (I18p . All the other details are as for 
Fig. [3 
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The image reconstruction technique we have implemented breaks down in the multiple 
scattering regime. In future work, we plan to explore corrections to the single scattering 
approximation. In addition, since electromagnetic waves in random media are, in general, 



^polarization within the framework of the 
8j. Finally, we note that a technique for 



polarized, we also intend to explore the effects o 
generalized vector radiative transport equation 
fluorescence imaging of mesoscopic objects has recently been reported [17|. It would thus 
be of interest to investigate the fluorescent analog of SSOT. 



Acknowledgment 

We are grateful to Guillaume Bal and Alexander Katsevich for valuable discussions. This 
work was supported by the NSF under the grant DMS-0554100 and by the NIH under the 
grant R01EB004832. 



APPENDIX: DERIVATION OF EQ. (1151) 

Substitution of (jSJ) into ffT2"j) with r = r 2 and s = s 2 results in the following expression 
for the single-scattered intensity: 



/ s (ri,si;r 2 ,s 2 ) = I Q A(s 2 ,Si) J ^t a (r)^(r 3 , r)g(r, n) 

x<5(u(r 2 - r) - s 2 )<5(u(r - n) - Si)d 3 r , (A.l) 

where we have introduced the notation u(r) = r/r. In the following analysis, the manip- 
ulation of delta functions is done in accordance with the theory of generalized functions 




We now make the change of variables r = r! + R, R = RR, and d 3 r = d 3 R = R 2 dRd 2 R. 
The integral over d 2 R is immediately evaluated and (lA.ip becomes 



J s (ri,si;r 2 ,s 2 ) = J ^4(s 2 ,Si) J g(r 2 , r x + Rs^g^t + Rs x , n) 

x/i s (n + i2si)5(u(r 2 i - Rsi) - s 2 )R 2 dR . (A.2) 



We then write the remaining delta-function as 
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5(u — s 2 ) = S((pu — ips 2 )S (cos 9^ — cos 9g 2 ) . (A. 3) 

Here u = r 2 \ — Rsi and 9 and <p are polar angles of the respective unit vectors. It is 
convenient to work in a reference frame whose z-axis coincides with the source-detector line. 
We then find that <p^ = <pg 1 ± n. Consequently, 

%u - <^§ 2 ) = 5 (1^ — ^s 2 1 — ?r) • (A.4) 

We next write 



5 (cos 9{, — cos 6*a 



(A.5) 



where 



f(R) 



r 2 i — Rcos9i 



COS^o . 



(A.6) 



y/r^ -2r 21 Rcos9 1 + R 2 
It can be verified that if 9\ + 62 > vr, the equation f(R) = has no positive roots. In the 

opposite limit, however, there is one positive root R = L\. Note that the lengths L\ and L 2 

are defined by ( floT) and illustrated in Fig. [TJ We thus have 



R I 5(f(R)) = @(^-9 1 -9 2 )Li 



2 5(R-L A 



Computation of the above derivative is straightforward and yields 



(A.7) 



\f{Li) 

Collecting everything, we arrive at 



Li 



r 21 L 2 sin 2 (6'i 



(A. 



R 2 S(u(r 21 - Rs x ) - s 2 ] 



S (|^si — VS2I — ®( n O2) 



r 2X L x L 2 b(R-L x \ 
sm 2 (9 x + 9 2 ) 



We then recall that ri + L±s = R21, L\ = |R 2 i — ri|, L 2 = |r 2 — R21I and obtain 



(A.9) 



g(r 2 ,r 1 + L^] 



1 



-2 exp 



l 2 



1 

g(rx + LiSi,ri) = exp 
^1 



/it (R 2 i + s 2 



A*t( r i + si 



(A.10) 
(A.11) 
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Finally, using (fl"5l) . we arrive at the result (fT3|) . 
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